Anàlisi de Components Principals (PCA)

1 Dades multivariants

Entenem per un conjunt multivariant de dades aquell en que s’han mesurat \(p\) variables \(X_1,...,X_p\) sobre un conjunt de \(n\) individus.

Els individus poden ser considerats com a punts d’un espai p-dimensional.

Les dades es poden representar en forma matricial

\[ \mathbf{X} = \left( \begin{array}{ccccc} x_{11} & x_{12} & x_{13} & \cdots & x_{1p} \\ x_{21} & x_{22} & x_{23} & \cdots & x_{2p} \\ \cdots & \cdots & \cdots & \cdots & \cdots \\ x_{n1} & x_{n2} & x_{n3} & \cdots & x_{np} \\ \end{array} \right)\]

on \(x_{ij}\) és el valor de la variable \(X_j\) sobre l’individu \(i\).

1.1 Introducció al PCA

• L’anàlisi de components principals (PCA) és una de les tècniques estadístiques de reducció de la dimensió més utilitzades.

• Destaca per la simplicitat del seu càlcul i per la intuïtiva interpretabilitat dels resultats obtinguts.

• Podem dir que l’objectiu principal és descobrir si hi han estructures amagades en un conjunt de dades, i en cercar aquest objectiu ens permet

– Identificar si diferents variables treballen juntes per crear propietats del conjunt de dades.

– Reduir la dimensió del conjunt de dades simplificant la descripció de les dades.

– Disminuir la redundància i el soroll en les dades.

– Preparar les dades per a posteriors anàlisis.

1.1.1 Principis bàsics del PCA

• Una forta correlació entre variables és un indicador d’una alta redundància.

• Els aspectes més rellevants del conjunt de dades són aquells que recullen el màxim de variabilitat possible.

1.2 PCA?

Una forma diferent de mirar les dades tractant d’esbrinar les principals característiques que potser queden amagades al primer cop de vista.

Quina imatge és més informativa?

El PCA va ser desenvolupat inicialment per Karl Pearson a finals del segle XIX i posteriorment estudiat per Hotelling als anys 30 del segle XX. Però no va poder popularitzar-se i arribar al nivell d’utilització actual fins l’aparició dels ordinadors.

1.3 Com funciona un PCA?

• L’objectiu és transformar el conjunt original de variables en un altre conjunt de noves variables incorrelacionades entre sí anomenat conjunt de components principals.

\[ X_1 , ..., X_p \Longrightarrow Y_1,...,Y_p \]

\[ Cov(Y_i,Y_j)=0 \,\,\,\,\, i \ne j\]

• Les noves variables són combinacions lineals de les anteriors i es van construint segons l’ordre en importància respecte la variabilitat total que recullen de la mostra.

\[ Y_i = \alpha_{i1} X_1 + \alpha_{i2} X_2 + \cdots + \alpha_{ip} X_p \,\,\,\,\, i=1,...,p\]

\[Var(Y_1) \ge Var(Y_2) \ge \cdots \ge Var(Y_p) \]

• Els coeficients s’assumeixen normalitzats per evitar incrementar la variabilitat artificialment.

\[\sum_{i=1}^p \alpha_i^2 = 1\]

• Idealment es tracta de poder seleccionar un conjunt de noves variables menor que el de variables originals i que estiguin incorrelacionades recollint la major part de la informació de les dades.

• La nova matriu de dades en els components principals és

\[ Y = \left( \begin{array}{ccccc} y_{11} & y_{12} & y_{13} & \cdots & y_{1p} \\ y_{21} & y_{22} & y_{23} & \cdots & y_{2p} \\ \cdots & \cdots & \cdots & \cdots & \cdots \\ y_{n1} & y_{n2} & y_{n3} & \cdots & y_{np} \\ \end{array} \right) = \left( \begin{array}{ccccc} x_{11} & x_{12} & x_{13} & \cdots & x_{1p} \\ x_{21} & x_{22} & x_{23} & \cdots & x_{2p} \\ \cdots & \cdots & \cdots & \cdots & \cdots \\ x_{n1} & x_{n2} & x_{n3} & \cdots & x_{np} \\ \end{array} \right) \cdot \left( \begin{array}{ccccc} \alpha_{11} & \alpha_{12} & \alpha_{13} & \cdots & \alpha_{1p} \\ \alpha_{21} & \alpha_{22} & \alpha_{23} & \cdots & \alpha_{2p} \\ \cdots & \cdots & \cdots & \cdots & \cdots \\ \alpha_{p1} & \alpha_{p2} & \alpha_{p3} & \cdots & \alpha_{pp} \\ \end{array} \right) = \mathbf{X} \cdot \mathbf{P} \]

on \(y_{ij} = \sum_{k=1}^p x_{ik} \alpha_{kj}\).

1.4 Exemple 1: Propietats químiques de diferents vins

En el fitxer wine.csv disposem de les dades d’una anàlisi química de vins provinents d’una mateixa regió de Itàlia però pertanyents a tres diferents varietats. L’analítica ha determinat la composició de 13 diferents constituents que són els següents:

1. Alcohol
2. Malic acid
3. Ash
4. Alcalinity of ash
5. Magnesium
6. Total phenols
7. Flavanoids
8. Nonflavanoid phenols
9. Proanthocyanins
10. Color intensity
11. Hue
12. OD280/OD315 of diluted wines
13. Proline

Fem primer una descriptiva de les dades

library(car)
wine<-read.table('../dades/wine.csv',header=T)
scatterplotMatrix(wine[,2:14],diagonal=list(method="boxplot"),smooth=F,groups=wine[,1])

round(cor(wine[,2:14]),2)

##             Alcohol Malic_ac   Ash Alcal_ash Magnesium phenols Flavanoids
## Alcohol        1.00     0.09  0.21     -0.31      0.27    0.29       0.24
## Malic_ac       0.09     1.00  0.16      0.29     -0.05   -0.34      -0.41
## Ash            0.21     0.16  1.00      0.44      0.29    0.13       0.12
## Alcal_ash     -0.31     0.29  0.44      1.00     -0.08   -0.32      -0.35
## Magnesium      0.27    -0.05  0.29     -0.08      1.00    0.21       0.20
## phenols        0.29    -0.34  0.13     -0.32      0.21    1.00       0.86
## Flavanoids     0.24    -0.41  0.12     -0.35      0.20    0.86       1.00
## Nonflav_phe   -0.16     0.29  0.19      0.36     -0.26   -0.45      -0.54
## Proanthocy     0.14    -0.22  0.01     -0.20      0.24    0.61       0.65
## Color_int      0.55     0.25  0.26      0.02      0.20   -0.06      -0.17
## Hue           -0.07    -0.56 -0.07     -0.27      0.06    0.43       0.54
## OD280.OD315    0.07    -0.37  0.00     -0.28      0.07    0.70       0.79
## Proline        0.64    -0.19  0.22     -0.44      0.39    0.50       0.49
##             Nonflav_phe Proanthocy Color_int   Hue OD280.OD315 Proline
## Alcohol           -0.16       0.14      0.55 -0.07        0.07    0.64
## Malic_ac           0.29      -0.22      0.25 -0.56       -0.37   -0.19
## Ash                0.19       0.01      0.26 -0.07        0.00    0.22
## Alcal_ash          0.36      -0.20      0.02 -0.27       -0.28   -0.44
## Magnesium         -0.26       0.24      0.20  0.06        0.07    0.39
## phenols           -0.45       0.61     -0.06  0.43        0.70    0.50
## Flavanoids        -0.54       0.65     -0.17  0.54        0.79    0.49
## Nonflav_phe        1.00      -0.37      0.14 -0.26       -0.50   -0.31
## Proanthocy        -0.37       1.00     -0.03  0.30        0.52    0.33
## Color_int          0.14      -0.03      1.00 -0.52       -0.43    0.32
## Hue               -0.26       0.30     -0.52  1.00        0.57    0.24
## OD280.OD315       -0.50       0.52     -0.43  0.57        1.00    0.31
## Proline           -0.31       0.33      0.32  0.24        0.31    1.00

1.5 Resum dels pasos per realitzar un PCA

• Construir la matriu de dades de forma que sigui centrada (a cada variable se li resta la seva mitjana) i si les escales de les variables són molt diferents entre elles (variàncies molt diferents) també és convenient escalar totes les variables dividint per la seva desviació típica per tal que totes tinguin variància 1.

• Si no escalem i una de les variables absorbeix molta variància tindrà molt pes en les components principals i sabem que una variància pot fer-se molt gran de forma artificial canviat les unitats de la variable.

• Calcular \(A=\mathbf{X}^T \mathbf{X}\)

– Si les variables són centrades la matriu \(A= \frac{1}{n-1} \mathbf{X}^T \mathbf{X}\) representa la matriu de covariàncies entre les variables.

– Si les variables també estan escalades, \(A= \frac{1}{n-1} \mathbf{X}^T \mathbf{X}\) representa la matriu de correlacions.

Per tant podem parlar de PCA basant en la matriu de covariàncies o PCA basat en la matriu de correlacions.

• Diagonalitzar la matriu \(A\).

– Els vectors propis són els coeficients dels components principals i els valors propis dividits per \(n-1\) ens donen les variàncies de cada component.

Podeu trobar més informació en les següents fons

• PCA en la Wikipedia

1.6 Propietats de les components principals

Sigui \(X\) una matriu de dades p-dimensional i siguin \(Y_1,...,Y_p\) els components principals

\[ Y_i = \alpha_{i1} X_1 + \alpha_{i2} X_2 + \cdots + \alpha_{ip} X_p \,\,\,\,\, i=1,...,p\]

on

\[\sum_{i=1}^p \alpha_i^2 = 1\]

Llavors

1. \(E(Y_i)=0\)

2. \(Cov(Y_i,Y_j)=0\)

3. \(Var(Y_1) \ge Var(Y_2) \ge ... \ge Var(Y_p)\)

4. Cap combinació lineal amb coeficients normalitzats \((\sum_{i=1}^p \alpha_i^2 = 1)\) de les variables \(X\) té una variància més gran que la de la primera component.

5. Donada una combinació lineal de coeficients normalitzats de les variables \(X\) incorrelacionada amb les \(k\) primeres components principals, la seva variància és màxima si es tracta de la component principal \(k+1\).

1.7 PCA amb R

• R porta incorporades en la distribució base dues funcions per realitzar un PCA: prcomp() i princomp().

• Les dues funcions són similars però prcomp() utilitza un algorisme una mica diferent (descomposició en valors singulars SVD) que ofereix una millor precisió numèrica i per tant és la funció preferida. princomp() no funciona si el número d’observacions (individus) és menor que el número de variables.

1.7.1 Exemple 1 amb prcomp()

wine.pca<-prcomp(wine[,2:14],scale=T) # per defecte scale=F,center=T
summary(wine.pca) # Importància de les components

## Importance of components:
##                          PC1    PC2    PC3     PC4     PC5     PC6     PC7
## Standard deviation     2.169 1.5802 1.2025 0.95863 0.92370 0.80103 0.74231
## Proportion of Variance 0.362 0.1921 0.1112 0.07069 0.06563 0.04936 0.04239
## Cumulative Proportion  0.362 0.5541 0.6653 0.73599 0.80162 0.85098 0.89337
##                            PC8     PC9   PC10    PC11    PC12    PC13
## Standard deviation     0.59034 0.53748 0.5009 0.47517 0.41082 0.32152
## Proportion of Variance 0.02681 0.02222 0.0193 0.01737 0.01298 0.00795
## Cumulative Proportion  0.92018 0.94240 0.9617 0.97907 0.99205 1.00000

plot(wine.pca)

wine.pca$rotation # Coeficients de les components

##                      PC1          PC2         PC3         PC4         PC5
## Alcohol     -0.144329395 -0.483651548 -0.20738262 -0.01785630  0.26566365
## Malic_ac     0.245187580 -0.224930935  0.08901289  0.53689028 -0.03521363
## Ash          0.002051061 -0.316068814  0.62622390 -0.21417556  0.14302547
## Alcal_ash    0.239320405  0.010590502  0.61208035  0.06085941 -0.06610294
## Magnesium   -0.141992042 -0.299634003  0.13075693 -0.35179658 -0.72704851
## phenols     -0.394660845 -0.065039512  0.14617896  0.19806835  0.14931841
## Flavanoids  -0.422934297  0.003359812  0.15068190  0.15229479  0.10902584
## Nonflav_phe  0.298533103 -0.028779488  0.17036816 -0.20330102  0.50070298
## Proanthocy  -0.313429488 -0.039301722  0.14945431  0.39905653 -0.13685982
## Color_int    0.088616705 -0.529995672 -0.13730621  0.06592568  0.07643678
## Hue         -0.296714564  0.279235148  0.08522192 -0.42777141  0.17361452
## OD280.OD315 -0.376167411  0.164496193  0.16600459  0.18412074  0.10116099
## Proline     -0.286752227 -0.364902832 -0.12674592 -0.23207086  0.15786880
##                     PC6         PC7         PC8         PC9        PC10
## Alcohol     -0.21353865 -0.05639636 -0.39613926 -0.50861912 -0.21160473
## Malic_ac    -0.53681385  0.42052391 -0.06582674  0.07528304  0.30907994
## Ash         -0.15447466 -0.14917061  0.17026002  0.30769445  0.02712539
## Alcal_ash    0.10082451 -0.28696914 -0.42797018 -0.20044931 -0.05279942
## Magnesium   -0.03814394  0.32288330  0.15636143 -0.27140257 -0.06787022
## phenols      0.08412230 -0.02792498  0.40593409 -0.28603452  0.32013135
## Flavanoids   0.01892002 -0.06068521  0.18724536 -0.04957849  0.16315051
## Nonflav_phe  0.25859401  0.59544729  0.23328465 -0.19550132 -0.21553507
## Proanthocy   0.53379539  0.37213935 -0.36822675  0.20914487 -0.13418390
## Color_int    0.41864414 -0.22771214  0.03379692 -0.05621752  0.29077518
## Hue         -0.10598274  0.23207564 -0.43662362 -0.08582839  0.52239889
## OD280.OD315 -0.26585107 -0.04476370  0.07810789 -0.13722690 -0.52370587
## Proline     -0.11972557  0.07680450 -0.12002267  0.57578611 -0.16211600
##                    PC11        PC12        PC13
## Alcohol      0.22591696  0.26628645 -0.01496997
## Malic_ac    -0.07648554 -0.12169604 -0.02596375
## Ash          0.49869142  0.04962237  0.14121803
## Alcal_ash   -0.47931378  0.05574287 -0.09168285
## Magnesium   -0.07128891 -0.06222011 -0.05677422
## phenols     -0.30434119  0.30388245  0.46390791
## Flavanoids   0.02569409  0.04289883 -0.83225706
## Nonflav_phe -0.11689586 -0.04235219 -0.11403985
## Proanthocy   0.23736257  0.09555303  0.11691707
## Color_int   -0.03183880 -0.60422163  0.01199280
## Hue          0.04821201 -0.25921400  0.08988884
## OD280.OD315 -0.04642330 -0.60095872  0.15671813
## Proline     -0.53926983  0.07940162 -0.01444734

head(wine.pca$x) # Coordenades dels individus en les noves components

##            PC1        PC2        PC3        PC4        PC5        PC6
## [1,] -3.307421 -1.4394023 -0.1652728 -0.2150246 -0.6910933 -0.2232504
## [2,] -2.203250  0.3324551 -2.0207571 -0.2905387  0.2569299 -0.9245123
## [3,] -2.509661 -1.0282507  0.9800541  0.7228632  0.2503270  0.5477310
## [4,] -3.746497 -2.7486184 -0.1756962  0.5663856  0.3109644  0.1141091
## [5,] -1.006070 -0.8673840  2.0209873 -0.4086131 -0.2976180 -0.4053761
## [6,] -3.041674 -2.1164309 -0.6276254 -0.5141870  0.6302409  0.1230834
##              PC7         PC8         PC9        PC10       PC11          PC12
## [1,]  0.59474883  0.06495586 -0.63963836 -1.01808396  0.4502932 -0.5392891439
## [2,]  0.05362434  1.02153432  0.30797798 -0.15925214  0.1422560 -0.3871456499
## [3,]  0.42301218 -0.34324787  1.17452129 -0.11304198  0.2858665 -0.0005819316
## [4,] -0.38225899  0.64178311 -0.05239662 -0.23873915 -0.7574476  0.2413387757
## [5,]  0.44282531  0.41552831 -0.32589984  0.07814604  0.5244656  0.2160546934
## [6,]  0.40052393  0.39378261  0.15171810  0.10170891 -0.4044444  0.3783653606
##              PC13
## [1,]  0.066052305
## [2,] -0.003626273
## [3,] -0.021655423
## [4,]  0.368444194
## [5,]  0.079140320
## [6,] -0.144747017

plot(wine.pca$x[,1],wine.pca$x[,2],main='PCA wine data set',
     xlab='PC1 (Flavonoides)',ylab='PC2 (Intensitat)',type='n')
text(wine.pca$x[,1],wine.pca$x[,2],label=wine[,1],cex=0.5,col=wine[,1])

library(rgl)
plot3d(wine.pca$x[,1],wine.pca$x[,2],wine.pca$x[,2],xlab='PC1',ylab='PC2',zlab='PC3',type='n')
text3d(wine.pca$x[,1],wine.pca$x[,2],wine.pca$x[,3],text=wine[,1],cex=0.5,col=wine[,1])

1.8 Biplot

• Un biplot és una representació simultània dels individus i les variables del conjunt de dades.

• Es pot aconseguir gràcies a la matriu que ens dóna el coeficient de cada variable original en cadascuna de les components principals

• La representació de les variables ajuda a interpretar les components principals

biplot(wine.pca,cex=c(0.3,1),xlim=c(-0.2,0.2)) # Representació simultània de individus i variables

1.9 PCA amb el package FactoMineR

library(FactoMineR)
wine$varietat<-as.factor(wine$varietat)
wine.pcaF<-PCA(wine,graph = F,quali.sup=1)
summary(wine.pcaF,nbelements=Inf,nbind=10,ncp=5)

## 
## Call:
## PCA(X = wine, quali.sup = 1, graph = F) 
## 
## 
## Eigenvalues
##                        Dim.1   Dim.2   Dim.3   Dim.4   Dim.5   Dim.6   Dim.7
## Variance               4.706   2.497   1.446   0.919   0.853   0.642   0.551
## % of var.             36.199  19.207  11.124   7.069   6.563   4.936   4.239
## Cumulative % of var.  36.199  55.406  66.530  73.599  80.162  85.098  89.337
##                        Dim.8   Dim.9  Dim.10  Dim.11  Dim.12  Dim.13
## Variance               0.348   0.289   0.251   0.226   0.169   0.103
## % of var.              2.681   2.222   1.930   1.737   1.298   0.795
## Cumulative % of var.  92.018  94.240  96.170  97.907  99.205 100.000
## 
## Individuals (the 10 first)
##                 Dist    Dim.1    ctr   cos2    Dim.2    ctr   cos2    Dim.3
## 1           |  4.000 |  3.317  1.313  0.687 |  1.443  0.469  0.130 | -0.166
## 2           |  3.384 |  2.209  0.583  0.426 | -0.333  0.025  0.010 | -2.026
## 3           |  3.322 |  2.517  0.756  0.574 |  1.031  0.239  0.096 |  0.983
## 4           |  4.855 |  3.757  1.685  0.599 |  2.756  1.709  0.322 | -0.176
## 5           |  2.668 |  1.009  0.122  0.143 |  0.870  0.170  0.106 |  2.027
## 6           |  3.945 |  3.050  1.111  0.598 |  2.122  1.013  0.289 | -0.629
## 7           |  3.384 |  2.449  0.716  0.524 |  1.175  0.311  0.121 | -0.977
## 8           |  3.348 |  2.059  0.506  0.378 |  1.609  0.582  0.231 |  0.146
## 9           |  3.528 |  2.511  0.753  0.507 |  0.918  0.190  0.068 | -1.771
## 10          |  3.274 |  2.754  0.905  0.707 |  0.789  0.140  0.058 | -0.984
##                ctr   cos2    Dim.4    ctr   cos2    Dim.5    ctr   cos2  
## 1            0.011  0.002 | -0.216  0.028  0.003 | -0.693  0.316  0.030 |
## 2            1.595  0.359 | -0.291  0.052  0.007 |  0.258  0.044  0.006 |
## 3            0.375  0.088 |  0.725  0.321  0.048 |  0.251  0.041  0.006 |
## 4            0.012  0.001 |  0.568  0.197  0.014 |  0.312  0.064  0.004 |
## 5            1.596  0.577 | -0.410  0.103  0.024 | -0.298  0.059  0.013 |
## 6            0.154  0.025 | -0.516  0.163  0.017 |  0.632  0.263  0.026 |
## 7            0.371  0.083 | -0.066  0.003  0.000 |  1.028  0.696  0.092 |
## 8            0.008  0.002 | -1.193  0.870  0.127 | -0.077  0.004  0.001 |
## 9            1.218  0.252 |  0.056  0.002  0.000 |  0.892  0.524  0.064 |
## 10           0.376  0.090 |  0.349  0.075  0.011 |  0.469  0.145  0.020 |
## 
## Variables
##                Dim.1    ctr   cos2    Dim.2    ctr   cos2    Dim.3    ctr
## Alcohol     |  0.313  2.083  0.098 |  0.764 23.392  0.584 | -0.249  4.301
## Malic_ac    | -0.532  6.012  0.283 |  0.355  5.059  0.126 |  0.107  0.792
## Ash         | -0.004  0.000  0.000 |  0.499  9.990  0.249 |  0.753 39.216
## Alcal_ash   | -0.519  5.727  0.270 | -0.017  0.011  0.000 |  0.736 37.464
## Magnesium   |  0.308  2.016  0.095 |  0.473  8.978  0.224 |  0.157  1.710
## phenols     |  0.856 15.576  0.733 |  0.103  0.423  0.011 |  0.176  2.137
## Flavanoids  |  0.917 17.887  0.842 | -0.005  0.001  0.000 |  0.181  2.271
## Nonflav_phe | -0.648  8.912  0.419 |  0.045  0.083  0.002 |  0.205  2.903
## Proanthocy  |  0.680  9.824  0.462 |  0.062  0.154  0.004 |  0.180  2.234
## Color_int   | -0.192  0.785  0.037 |  0.837 28.090  0.701 | -0.165  1.885
## Hue         |  0.644  8.804  0.414 | -0.441  7.797  0.195 |  0.102  0.726
## OD280.OD315 |  0.816 14.150  0.666 | -0.260  2.706  0.068 |  0.200  2.756
## Proline     |  0.622  8.223  0.387 |  0.577 13.315  0.332 | -0.152  1.606
##               cos2    Dim.4    ctr   cos2    Dim.5    ctr   cos2  
## Alcohol      0.062 | -0.017  0.032  0.000 |  0.245  7.058  0.060 |
## Malic_ac     0.011 |  0.515 28.825  0.265 | -0.033  0.124  0.001 |
## Ash          0.567 | -0.205  4.587  0.042 |  0.132  2.046  0.017 |
## Alcal_ash    0.542 |  0.058  0.370  0.003 | -0.061  0.437  0.004 |
## Magnesium    0.025 | -0.337 12.376  0.114 | -0.672 52.860  0.451 |
## phenols      0.031 |  0.190  3.923  0.036 |  0.138  2.230  0.019 |
## Flavanoids   0.033 |  0.146  2.319  0.021 |  0.101  1.189  0.010 |
## Nonflav_phe  0.042 | -0.195  4.133  0.038 |  0.463 25.070  0.214 |
## Proanthocy   0.032 |  0.383 15.925  0.146 | -0.126  1.873  0.016 |
## Color_int    0.027 |  0.063  0.435  0.004 |  0.071  0.584  0.005 |
## Hue          0.011 | -0.410 18.299  0.168 |  0.160  3.014  0.026 |
## OD280.OD315  0.040 |  0.177  3.390  0.031 |  0.093  1.023  0.009 |
## Proline      0.023 | -0.222  5.386  0.049 |  0.146  2.492  0.021 |
## 
## Supplementary categories
##                  Dist     Dim.1    cos2  v.test     Dim.2    cos2  v.test  
## varietat_1  |   2.515 |   2.283   0.824   9.858 |   0.968   0.148   5.738 |
## varietat_2  |   1.679 |  -0.039   0.001  -0.195 |  -1.643   0.959 -11.272 |
## varietat_3  |   3.030 |  -2.748   0.823 -10.242 |   1.241   0.168   6.351 |
##               Dim.3    cos2  v.test     Dim.4    cos2  v.test     Dim.5    cos2
## varietat_1   -0.160   0.004  -1.243 |  -0.112   0.002  -1.098 |   0.243   0.009
## varietat_2    0.262   0.024   2.358 |   0.051   0.001   0.576 |  -0.122   0.005
## varietat_3   -0.191   0.004  -1.283 |   0.063   0.000   0.529 |  -0.118   0.002
##              v.test  
## varietat_1    2.460 |
## varietat_2   -1.432 |
## varietat_3   -1.030 |

Per obtenir els coeficients dels components (loadings)

sweep(wine.pcaF$var$coord,2,sqrt(wine.pcaF$eig[1:ncol(wine.pcaF$var$coord),1]),FUN="/")

##                    Dim.1        Dim.2       Dim.3       Dim.4       Dim.5
## Alcohol      0.144329395  0.483651548 -0.20738262 -0.01785630  0.26566365
## Malic_ac    -0.245187580  0.224930935  0.08901289  0.53689028 -0.03521363
## Ash         -0.002051061  0.316068814  0.62622390 -0.21417556  0.14302547
## Alcal_ash   -0.239320405 -0.010590502  0.61208035  0.06085941 -0.06610294
## Magnesium    0.141992042  0.299634003  0.13075693 -0.35179658 -0.72704851
## phenols      0.394660845  0.065039512  0.14617896  0.19806835  0.14931841
## Flavanoids   0.422934297 -0.003359812  0.15068190  0.15229479  0.10902584
## Nonflav_phe -0.298533103  0.028779488  0.17036816 -0.20330102  0.50070298
## Proanthocy   0.313429488  0.039301722  0.14945431  0.39905653 -0.13685982
## Color_int   -0.088616705  0.529995672 -0.13730621  0.06592568  0.07643678
## Hue          0.296714564 -0.279235148  0.08522192 -0.42777141  0.17361452
## OD280.OD315  0.376167411 -0.164496193  0.16600459  0.18412074  0.10116099
## Proline      0.286752227  0.364902832 -0.12674592 -0.23207086  0.15786880

plot(wine.pcaF,choix='ind',habillage='ind',col.hab=wine$varietat,axes=c(1,2))

plot(wine.pcaF,choix='var')

plotellipses(wine.pcaF,main='Confidence for the means of the categories')

plotellipses(wine.pcaF,means=F)

Per aconseguir un biplot

#install.packages('devtools')
library(devtools)
install_github("fawda123/ggord")
library(ggord)
ggord(wine.pcaF,wine$varietat,poly=F,size=2,ext=1.5,repel=T)

Un nou package que realitza un estudi previ amb molta informació, però que ha de ser contrastat amb l’anàlisi acurada dels resultats, és el package FactoInvestigate, automàticament crea un informe, per defecte en format html, sobre els resultats del PCA.

library(FactoInvestigate)
Investigate(wine.pcaF)

1.10 Exemple 2: Composició d’aminoàcids als teixits animals

Les dades corresponen als percentatges dels principals aminoàcids en el teixit muscular de diferents animals i en sis òrgans interns dels bous (beef).

amino<-read.table('../dades/Aminoacids_Composition.csv',sep=';',dec=',',header=T)
head(amino)

Animal_tissue	Arginine	Histidine	Lysine	Phenylalanine	Tyrosisne	Tryptophane	Serine	Threonine	Cystine_G	Cystine__S	Methionine_B	Methionine_S	Total_S	Organic_S
Beef	6.91	2.25	8.11	4.92	4.30	1.35	5.43	4.57	1.29	0.94	3.17	3.11	1.09	1.07
Veal	7.54	2.39	9.62	4.41	4.86	1.39	6.06	5.10	1.34	0.99	3.28	3.62	1.14	1.12
Lamb	7.59	2.37	8.68	4.47	4.89	1.44	6.29	5.28	1.42	0.99	3.06	3.28	1.15	1.13
Pork	6.62	2.16	8.65	3.95	4.41	1.31	4.57	4.61	1.14	0.99	3.42	3.22	1.03	1.02
Chickenlight	6.91	2.34	8.44	3.85	4.23	1.30	4.72	4.66	1.26	0.82	3.28	3.35	1.06	1.04
Chicken_dark	7.08	2.28	8.42	4.61	4.34	1.22	5.45	4.58	1.28	0.88	2.95	3.64	1.05	1.04

amino$Animal_tissue<-as.factor(amino$Animal_tissue)

amino.pca<-prcomp(amino[,2:15],scale=T)
summary(amino.pca)

## Importance of components:
##                           PC1    PC2    PC3     PC4     PC5     PC6     PC7
## Standard deviation     2.1990 2.1129 1.3930 0.87138 0.84740 0.68947 0.56761
## Proportion of Variance 0.3454 0.3189 0.1386 0.05424 0.05129 0.03396 0.02301
## Cumulative Proportion  0.3454 0.6643 0.8029 0.85712 0.90841 0.94236 0.96538
##                            PC8     PC9    PC10    PC11    PC12    PC13    PC14
## Standard deviation     0.41748 0.35314 0.30142 0.21311 0.17748 0.12852 0.03797
## Proportion of Variance 0.01245 0.00891 0.00649 0.00324 0.00225 0.00118 0.00010
## Cumulative Proportion  0.97783 0.98673 0.99322 0.99647 0.99872 0.99990 1.00000

plot(amino.pca)

amino.pca$rotation # Coeficients de les components

##                       PC1         PC2         PC3         PC4         PC5
## Arginine      -0.16729959  0.02031698 -0.46472267  0.55297031 -0.40830017
## Histidine     -0.30557250  0.13832749 -0.34010241 -0.34014596  0.31818090
## Lysine        -0.30579819 -0.31075579 -0.14538434  0.13089593  0.10854675
## Phenylalanine -0.17521244  0.36549740  0.21946471  0.07304489 -0.02334372
## Tyrosisne     -0.37961965  0.18331259  0.02973021  0.05971914 -0.05581050
## Tryptophane   -0.24229576  0.34675626  0.11533755 -0.07108400 -0.03625000
## Serine         0.16323709  0.33669835 -0.17974678 -0.14877137 -0.56398843
## Threonine     -0.31462716  0.21759569 -0.32714111 -0.09198794 -0.10191760
## Cystine_G     -0.05062144  0.43168253  0.01013767 -0.26518119  0.19406365
## Cystine__S     0.01663731  0.35531419  0.10953609  0.53624996  0.35895804
## Methionine_B  -0.36634438 -0.19425405 -0.14192627  0.11451349  0.28567177
## Methionine_S  -0.30506427 -0.27159985 -0.01298282 -0.34921605 -0.21976985
## Total_S       -0.29125771 -0.03230464  0.50708714 -0.02150397 -0.24457554
## Organic_S     -0.33320661 -0.10293426  0.39496212  0.15306390 -0.17489889
##                       PC6          PC7          PC8         PC9        PC10
## Arginine       0.39628966 -0.152848342 -0.146754887  0.11096139 -0.17292465
## Histidine      0.25101870 -0.340779026 -0.094611835 -0.18589570  0.10886175
## Lysine        -0.03942168 -0.141648386  0.145218889  0.21874383  0.45961881
## Phenylalanine  0.23593971  0.449828928 -0.464328716  0.09716822  0.54186244
## Tyrosisne     -0.43805108 -0.008736903 -0.136895831  0.54245348 -0.30351483
## Tryptophane    0.34695983  0.366679554  0.359068474 -0.13704005 -0.41560580
## Serine        -0.27322158 -0.063339105  0.154306598 -0.22207710  0.23868849
## Threonine     -0.35398683  0.114177918  0.244174693 -0.15511137  0.18614345
## Cystine_G      0.16478767 -0.390453857  0.005101925  0.25358050 -0.05486455
## Cystine__S    -0.34829358 -0.208262068 -0.203370841 -0.39643059 -0.08539546
## Methionine_B  -0.10466689  0.375742552  0.235940329 -0.12950680 -0.02667420
## Methionine_S  -0.13292832  0.037185555 -0.597144990 -0.30865661 -0.26942363
## Total_S       -0.05648098 -0.246464710  0.168507470  0.15746068  0.05502402
## Organic_S      0.17383134 -0.295602190  0.140317513 -0.38875602  0.10139067
##                      PC11        PC12         PC13        PC14
## Arginine      -0.18042968 -0.07133822  0.030061426 -0.02091773
## Histidine     -0.08377629  0.29496254 -0.154260742  0.45007220
## Lysine         0.53299285 -0.40777219 -0.043562986  0.07617816
## Phenylalanine -0.04838190  0.04550188 -0.016428336 -0.01773564
## Tyrosisne      0.20395527  0.37572522 -0.184656018  0.02779791
## Tryptophane    0.36888119 -0.21290278  0.028583385  0.20950379
## Serine         0.07889758 -0.09322749 -0.511566531  0.04564824
## Threonine     -0.13802973  0.06158054  0.654571698 -0.14074669
## Cystine_G     -0.08890417 -0.34984117 -0.046948839 -0.56882135
## Cystine__S     0.05974975 -0.22446354 -0.005374916  0.14348648
## Methionine_B  -0.42126250 -0.09485963 -0.488802084 -0.25003096
## Methionine_S   0.08573607 -0.32041399  0.030335668 -0.09122078
## Total_S       -0.48550293 -0.27915532  0.063903483  0.39853820
## Organic_S      0.19762720  0.42215509 -0.043457606 -0.38795147

head(amino.pca$x) # Coordenades dels individus en les noves components

##               PC1        PC2        PC3        PC4        PC5         PC6
## [1,] -0.068020724 -0.3112500 -0.5433524  0.3072096  0.5370535  0.58415098
## [2,] -2.556325380 -0.3066863 -1.6126095  0.7105575 -0.6274207 -0.24271605
## [3,] -2.187933293  0.6408842 -1.4432913  0.7932065 -0.9209551 -0.09590444
## [4,]  0.272296537 -1.5327358 -1.0841421  0.3610637  1.6211071 -0.51155039
## [5,] -0.059486256 -1.6384165 -1.5180538 -0.3831723  0.9622341  0.50473727
## [6,]  0.009382133 -1.0229977 -1.4119322 -0.1933520  0.1115477  0.34866612
##              PC7          PC8         PC9         PC10        PC11        PC12
## [1,]  0.42984859  0.042729826 -0.06097912  0.574771343 -0.17376331  0.06311394
## [2,] -0.74503158 -0.008727959  0.04749187  0.095201640  0.27111662 -0.25103512
## [3,] -0.89087170  0.447963433  0.18723414  0.007258855 -0.01311947  0.22301517
## [4,]  0.79385438  0.350037839 -0.01027192 -0.290582805  0.30749002 -0.07186221
## [5,]  0.08217864  0.364275648 -0.04325920 -0.218347390 -0.17365500 -0.02135495
## [6,]  0.08092372 -0.957533023 -0.01345458  0.246912006  0.16627456 -0.17896378
##             PC13         PC14
## [1,] -0.17016655  0.018299851
## [2,] -0.17682681  0.060249494
## [3,]  0.20787711 -0.032654909
## [4,]  0.04165848  0.019623646
## [5,]  0.19139162  0.003231091
## [6,]  0.04873566 -0.090452445

plot(amino.pca$x[,1],amino.pca$x[,2],main='PCA amino acids data set',
     xlab='PC1',ylab='PC2',type='n')
text(amino.pca$x[,1],amino.pca$x[,2],label=amino[,1],cex=0.5)

biplot(amino.pca,xlim=c(-0.5,0.8)) # Representació simultània de individus i variables

1.11 Exemple 3: Contaminació a les ciutats d’Estats Units

Una anàlisi de components principals pot ser un pas previ per aplicar una altra tècnica estadística. Per exemple es molt freqüent aplicar un PCA a unes dades i posteriorment una regressió sobre les components obtingudes.

En el següent exemple utilitzarem unes dades incorporades al package HSAUR2 sobre contaminació en diverses ciutats americanes. La variable dependent d’interès és la concentració de \(SO_2\) i la resta són diverses variables predictores.

Imaginem que volem fer una regressió lineal múltiple

library(HSAUR2)
library(car)
data(USairpollution)
head(USairpollution)

	SO2	temp	manu	popul	wind	precip	predays
Albany	46	47.6	44	116	8.8	33.36	135
Albuquerque	11	56.8	46	244	8.9	7.77	58
Atlanta	24	61.5	368	497	9.1	48.34	115
Baltimore	47	55.0	625	905	9.6	41.31	111
Buffalo	11	47.1	391	463	12.4	36.11	166
Charleston	31	55.2	35	71	6.5	40.75	148

scatterplotMatrix(USairpollution,diag='boxplot')

air.lm<-lm(SO2~.,data=USairpollution)
summary(air.lm)

## 
## Call:
## lm(formula = SO2 ~ ., data = USairpollution)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -23.004  -8.542  -0.991   5.758  48.758 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 111.72848   47.31810   2.361 0.024087 *  
## temp         -1.26794    0.62118  -2.041 0.049056 *  
## manu          0.06492    0.01575   4.122 0.000228 ***
## popul        -0.03928    0.01513  -2.595 0.013846 *  
## wind         -3.18137    1.81502  -1.753 0.088650 .  
## precip        0.51236    0.36276   1.412 0.166918    
## predays      -0.05205    0.16201  -0.321 0.749972    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 14.64 on 34 degrees of freedom
## Multiple R-squared:  0.6695, Adjusted R-squared:  0.6112 
## F-statistic: 11.48 on 6 and 34 DF,  p-value: 5.419e-07

Però, compte!, hem de vigilar les premisses, entre elles, verifiquem si tenim problemes de multicol·linealitat

vif(air.lm)

##      temp      manu     popul      wind    precip   predays 
##  3.763996 14.703652 14.340833  1.255519  3.404921  3.443651

Què fem? Eliminem alguna de les variables? Tenim alguna altra alternativa?

1.11.1 Sí, la regressió en components principals

La regressió en components principals (PCR) combina l’anàlisi de components principals (PCA) amb la regressió lineal per manejar multicol·linealitat i alta dimensionalitat. Primer, el PCA transforma les variables predictives correlacionades en components principals no correlacionats, ordenats per variància explicada. Aquests components s’utilitzen com a predictors en un model de regressió per predir la variable resposta. Aquest mètode permet redueix la dimensionalitat, millora l’estabilitat i evita el sobreajust, ideal per dades òmiques o multivariantes.

air.pca<-prcomp(USairpollution[,-1],scale=T)
air.pca$rotation

##                 PC1        PC2         PC3         PC4         PC5         PC6
## temp     0.32964613 -0.1275974  0.67168611  0.30645728  0.55805638  0.13618780
## manu    -0.61154243 -0.1680577  0.27288633 -0.13684076 -0.10204211  0.70297051
## popul   -0.57782195 -0.2224533  0.35037413 -0.07248126  0.07806551 -0.69464131
## wind    -0.35383877  0.1307915 -0.29725334  0.86942583  0.11326688  0.02452501
## precip   0.04080701  0.6228578  0.50456294  0.17114826 -0.56818342 -0.06062222
## predays -0.23791593  0.7077653 -0.09308852 -0.31130693  0.58000387  0.02196062

airnew<-data.frame(USairpollution$SO2,air.pca$x)
head(airnew)

	USairpollution.SO2	PC1	PC2	PC3	PC4	PC5	PC6
Albany	46	0.5323325	0.7823504	-1.3458855	-0.8719022	-0.0456667	-0.0618070
Albuquerque	11	1.3997048	-2.8307246	-1.2597507	0.0948076	0.2397712	0.0302114
Atlanta	24	0.5916448	0.5800300	0.9831972	0.2263505	-0.1166496	0.0587422
Baltimore	47	-0.5031297	0.0284011	0.3591380	0.0863224	-0.3186253	-0.1910289
Buffalo	11	-1.3736443	1.8572290	-1.7543973	0.8459397	0.7306775	0.0187745
Charleston	31	1.4122097	1.1957293	0.0784687	-1.9867817	0.2821117	0.0573742

Ara fem la regressió sobre les components principals

colnames(airnew)[1]<-'SO2'
airnew.lm<-lm(SO2~.,data=airnew)

Evidentment aquí no hi haurà cap problema de multicol·linealitat, recordem que totes les components principals són independents, veiem-lo

vif(airnew.lm)

## PC1 PC2 PC3 PC4 PC5 PC6 
##   1   1   1   1   1   1

summary(airnew.lm)

## 
## Call:
## lm(formula = SO2 ~ ., data = airnew)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -23.004  -8.542  -0.991   5.758  48.758 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  30.0488     2.2858  13.146 6.91e-15 ***
## PC1         -10.0655     1.5616  -6.446 2.28e-07 ***
## PC2           2.2674     1.8895   1.200  0.23845    
## PC3           0.3797     1.9596   0.194  0.84752    
## PC4          -8.6553     2.6541  -3.261  0.00253 ** 
## PC5         -15.3644     6.8369  -2.247  0.03122 *  
## PC6          39.7591    12.4686   3.189  0.00306 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 14.64 on 34 degrees of freedom
## Multiple R-squared:  0.6695, Adjusted R-squared:  0.6112 
## F-statistic: 11.48 on 6 and 34 DF,  p-value: 5.419e-07

Interpreteu la variable PCA6?

Podem inclús simplificar el nostre model

step(airnew.lm)

## Start:  AIC=226.37
## SO2 ~ PC1 + PC2 + PC3 + PC4 + PC5 + PC6
## 
##        Df Sum of Sq     RSS    AIC
## - PC3   1       8.0  7291.3 224.42
## - PC2   1     308.4  7591.7 226.07
## <none>               7283.3 226.37
## - PC5   1    1081.8  8365.1 230.05
## - PC6   1    2178.1  9461.4 235.10
## - PC4   1    2278.1  9561.3 235.53
## - PC1   1    8900.1 16183.4 257.11
## 
## Step:  AIC=224.42
## SO2 ~ PC1 + PC2 + PC4 + PC5 + PC6
## 
##        Df Sum of Sq     RSS    AIC
## - PC2   1     308.4  7599.8 224.11
## <none>               7291.3 224.42
## - PC5   1    1081.8  8373.2 228.09
## - PC6   1    2178.1  9469.5 233.13
## - PC4   1    2278.1  9569.4 233.56
## - PC1   1    8900.1 16191.4 255.12
## 
## Step:  AIC=224.11
## SO2 ~ PC1 + PC4 + PC5 + PC6
## 
##        Df Sum of Sq     RSS    AIC
## <none>               7599.8 224.11
## - PC5   1    1081.8  8681.6 227.57
## - PC6   1    2178.1  9777.9 232.45
## - PC4   1    2278.1  9877.8 232.86
## - PC1   1    8900.1 16499.9 253.90

## 
## Call:
## lm(formula = SO2 ~ PC1 + PC4 + PC5 + PC6, data = airnew)
## 
## Coefficients:
## (Intercept)          PC1          PC4          PC5          PC6  
##      30.049      -10.065       -8.655      -15.364       39.759

1.11.1.1 Avantatges i limitacions de fer regressió amb components principals

Avantatges

Redueix la dimensionalitat projectant dades en components no correlacionats, evitant multicol·linealitat i estabilitzant estimacions de coeficients. Millora la predicció en dades sorolloses seleccionant components amb major variància explicada, reduint sobreajust. És computacionalment eficient i escalable per grans datasets com dades òmiques.

Limitacions

Els components principals poden no capturar relacions no lineals ni interpretar-se fàcilment, perdent significat biològic o físic directe. La selecció subjectiva de components pot introduir biaixos, i assumeix linealitat entre components i resposta. Menys flexible que mètodes regulars com ridge o Lasso per control precís de penalització.

1.12 Exemple 4: Marques en l’heptatló olímpic

Les dades corresponen als resultats del heptatló olímpic als jocs olímpics de Seül al 1988. La última columna és la suma de les puntuacions assignades a cada prova per les fórmules estàndard establertes per World Athletics des de 1985, específiques per a cada disciplina. i per tant la classificació final de la prova.

Les fórmules són:

Fórmula general

Per esdeveniments de velocitat (100 m tanques, 200 m, 800 m): \(P = INT(A - B \cdot t^C)\), on \(t\) és el temps en segons, \(INT\) arrodoneix a l’enter més proper, i \(A, B, C\) són constants per prova. Per salts i llançaments (salt alçada/llargada, pes, javelina): \(P = INT(A \cdot (d - B)^C)\), on \(d\) és la distància en metres (salts en cm per alçada).

Constants per disciplina femenina

• 100 m tanques: \(A=9.230, B=26.7, C=1.835\)
  
• Salt alçada: \(A=1.845, B=0.688, C=1.6\)
  
• Llançament pes: \(A=56.021, B=4.99087, C=1.81\)
  
• 200 m: \(A=4.99087, B=42.5, C=1.81\)
  
• Salt llargada: \(A=0.188807, B=0.296, C=1.41\)
  
• Llançament javelina: \(A=15.9803, B=3.8, C=1.04\)
  
• 800 m: \(A=361.539, B=0.11193, C=1.88\).

library(HSAUR2)
data(heptathlon)
head(heptathlon)

	hurdles	highjump	shot	run200m	longjump	javelin	run800m	score
Joyner-Kersee (USA)	12.69	1.86	15.80	22.56	7.27	45.66	128.51	7291
John (GDR)	12.85	1.80	16.23	23.65	6.71	42.56	126.12	6897
Behmer (GDR)	13.20	1.83	14.20	23.10	6.68	44.54	124.20	6858
Sablovskaite (URS)	13.61	1.80	15.23	23.92	6.25	42.78	132.24	6540
Choubenkova (URS)	13.51	1.74	14.76	23.93	6.32	47.46	127.90	6540
Schulz (GDR)	13.75	1.83	13.50	24.65	6.33	42.82	125.79	6411

hept.pca<-prcomp(heptathlon[,-8],scale=T)
hept.pca$rotation

##                 PC1         PC2         PC3         PC4         PC5         PC6
## hurdles   0.4528710 -0.15792058  0.04514996 -0.02653873 -0.09494792 -0.78334101
## highjump -0.3771992  0.24807386 -0.36777902  0.67999172 -0.01879888 -0.09939981
## shot     -0.3630725 -0.28940743  0.67618919  0.12431725 -0.51165201  0.05085983
## run200m   0.4078950  0.26038545 -0.08359211  0.36106580 -0.64983404  0.02495639
## longjump -0.4562318  0.05587394  0.13931653  0.11129249  0.18429810 -0.59020972
## javelin  -0.0754090 -0.84169212 -0.47156016  0.12079924 -0.13510669  0.02724076
## run800m   0.3749594 -0.22448984  0.39585671  0.60341130  0.50432116  0.15555520
##                  PC7
## hurdles  -0.38024707
## highjump -0.43393114
## shot     -0.21762491
## run200m   0.45338483
## longjump  0.61206388
## javelin   0.17294667
## run800m   0.09830963

summary(hept.pca)

## Importance of components:
##                           PC1    PC2     PC3     PC4     PC5     PC6    PC7
## Standard deviation     2.1119 1.0928 0.72181 0.67614 0.49524 0.27010 0.2214
## Proportion of Variance 0.6372 0.1706 0.07443 0.06531 0.03504 0.01042 0.0070
## Cumulative Proportion  0.6372 0.8078 0.88223 0.94754 0.98258 0.99300 1.0000

biplot(hept.pca,xlim=c(-0.8,0.8),cex=c(0.55,1))

cor(hept.pca$x[,1],heptathlon$score)

## [1] -0.9910978

plot(hept.pca$x[,1],heptathlon$score,type='p')
text(hept.pca$x[,1],heptathlon$score,label=row.names(heptathlon),cex=0.5,pos=4)

No l’hauríem fet tant malament si haguéssim donat les medalles amb un criteri estadístic, no?

1.13 Exemple 5: Classificació de marques de Whisky

Les dades corresponen a les puntuacions obtingudes per diferents marques de whisky en diferents sabors a partir de l’opinió d’un panell d’experts.

whisky<-read.table('../dades/whiskies.txt',header=T,sep=',',dec='.')
head(whisky)

RowID	Distillery	Body	Sweetness	Smoky	Medicinal	Tobacco	Honey	Spicy	Winey	Nutty	Malty	Fruity	Floral	Postcode	Latitude	Longitude
1	Aberfeldy	2	2	2	0	0	2	1	2	2	2	2	2	PH15 2EB	286580	749680
2	Aberlour	3	3	1	0	0	4	3	2	2	3	3	2	AB38 9PJ	326340	842570
3	AnCnoc	1	3	2	0	0	2	0	0	2	2	3	2	AB5 5LI	352960	839320
4	Ardbeg	4	1	4	4	0	0	2	0	1	2	1	0	PA42 7EB	141560	646220
5	Ardmore	2	2	2	0	0	1	1	1	2	3	1	1	AB54 4NH	355350	829140
6	ArranIsleOf	2	3	1	1	0	1	1	1	0	1	1	2	KA27 8HJ	194050	649950

whisky.pca<-prcomp(whisky[,3:14])
summary(whisky.pca)

## Importance of components:
##                           PC1    PC2    PC3    PC4     PC5     PC6     PC7
## Standard deviation     1.5358 1.2270 0.8654 0.8039 0.75261 0.68513 0.63256
## Proportion of Variance 0.3011 0.1922 0.0956 0.0825 0.07231 0.05992 0.05108
## Cumulative Proportion  0.3011 0.4933 0.5889 0.6714 0.74370 0.80363 0.85471
##                            PC8     PC9    PC10    PC11    PC12
## Standard deviation     0.59943 0.52347 0.50049 0.42422 0.27266
## Proportion of Variance 0.04587 0.03498 0.03198 0.02297 0.00949
## Cumulative Proportion  0.90058 0.93556 0.96754 0.99051 1.00000

plot(whisky.pca)

whisky.pca$rotation[,1:5]

##                   PC1         PC2           PC3         PC4          PC5
## Body       0.36119005  0.49130643  0.0301178096  0.07460952 -0.227019231
## Sweetness -0.20298238  0.04659634 -0.2638792230  0.37063897 -0.009220720
## Smoky      0.47794419  0.06874216  0.2188100647 -0.08852143  0.201614612
## Medicinal  0.57527678 -0.16079484  0.0431598432 -0.08237288  0.033120119
## Tobacco    0.09173306 -0.02004776 -0.0006685359  0.03337611  0.008967789
## Honey     -0.22090804  0.41799491  0.1102471134 -0.03315990  0.596888644
## Spicy      0.05811101  0.17548310  0.6992443906  0.17163088  0.133792990
## Winey     -0.03745608  0.63964979 -0.2331295871  0.22573776 -0.110928262
## Nutty     -0.04766410  0.26036122 -0.1785529039 -0.85059012 -0.025288924
## Malty     -0.12781608  0.10296202  0.1084169899 -0.07177790  0.105401421
## Fruity    -0.20235755  0.12374977  0.4034627791 -0.09457148 -0.702770800
## Floral    -0.38394443 -0.13074914  0.3433008365 -0.14903423  0.120141367

biplot(whisky.pca,cex=c(0.1,1))

plot(whisky.pca$x[,1],whisky.pca$x[,2],col='white',
     xlab='C.P.1 (fumat)',ylab='C.P.2 (potent)')
text(whisky.pca$x[,1],whisky.pca$x[,2],label=whisky[,2],cex=0.5,col='blue')

1.14 Annex I: Demostració matemàtica de les components principals (Opcional)

1.14.1 Introducció

Donada una matriu de dades \(\mathbf{X}\), l’objectiu es trobar una nova matriu \(\mathbf{P}\) tal que la matriu de covariàncies de \(\mathbf{X}\mathbf{P}\) sigui diagonal i els termes de la diagonal estiguin en ordre descendent. Si ho aconseguim

• \(\mathbf{X}\mathbf{P}\) representarà el nou conjunt de dades

• Les noves variables seran combinacions lineals de les variables originals amb coeficients donats per les columnes de \(\mathbf{P}\)

• Si la matriu de covariàncies és diagonal la correlació entre cada parell de noves variables serà 0 i la variància serà el valor del terme de la diagonal.

• En el següent apartat es demostra que aquest objectiu s’assoleix diagonalitzant la matriu de covariàncies de les dades.

– Les columnes de P corresponen als vectors propis de la diagonalització

– Els valors propis ens proporcionen les variàncies de les noves variables, si els ordenem de més gran a més petit tindrem ordenades les components principals segons la variància explicada

1.14.2 Repàs de conceptes d’àlgebra lineal

• La matriu de covariàncies d’una matriu de dades centrades \(\mathbf{X}\) de dimensió \(n \times p\) és una matriu quadrada simètrica de dimensió \(p\) que ve donada per

\[ \frac{1}{n-1} \mathbf{X}^T \mathbf{X} \]

• Tota matriu simètrica \(\mathbf{A}\) d’ordre \(n\) té un conjunt de \(n\) vectors propis ortonormals \((v_1,...,v_n)\) amb valors propis associats \((\lambda_1,...,\lambda_n)\)

• Si \(\mathbf{A}\) és una matriu simètrica i \(\mathbf{E}\) és una matriu que té per columna i-èsima el i-èsim vector propi de \(\mathbf{A}\), llavors existeix una matriu diagonal \(\mathbf{D}\), tal que \(\mathbf{A}=\mathbf{E}\mathbf{D}\mathbf{E}^T\).

• Si la matriu \(\mathbf{E}\) és ortogonal (ortonormal) llavors \(\mathbf{E}^{-1}=\mathbf{E}^T \Rightarrow \mathbf{E} \cdot \mathbf{E}^T = \mathbf{E}^T \cdot \mathbf{E} = I\).

1.14.3 Procediment per obtenir P

\[ \begin{array}{llr} Cov(\mathbf{X}\mathbf{P}) & = \frac{1}{n-1} (\mathbf{X}\mathbf{P})^T(\mathbf{X}\mathbf{P}) & \\ & = \frac{1}{n-1} \mathbf{P}^T \mathbf{X}^T \mathbf{X} \mathbf{P} & \\ & = \frac{1}{n-1} \mathbf{P}^T \mathbf{A} \mathbf{P} & \,\,\,\,\,\,\,\,\,\, \mathbf{A}=\mathbf{X}^T \mathbf{X}\\ & = \frac{1}{n-1} \mathbf{P}^T \mathbf{E} \mathbf{D} \mathbf{E}^T \mathbf{P} & \,\,\,\,\,\,\,\,\,\, \mbox{veure apartat anterior}\\ & = \frac{1}{n-1} \mathbf{P}^T (\mathbf{P}\mathbf{D}\mathbf{P}^T )\mathbf{P} & \,\,\,\,\,\,\,\,\,\, \mathbf{P}=\mathbf{E}\\ & = \frac{1}{n-1} \mathbf{P}^{-1}\mathbf{P}\mathbf{D}\mathbf{P}^{-1}\mathbf{P} & \,\,\,\,\,\,\,\,\,\, \mathbf{P}^{-1}=\mathbf{P}^T\\ & = \frac{1}{n-1} \mathbf{D} & \end{array} \]

1.15 Annex II: Resolució amb les funcions de R bàsiques (Opcional)

1.15.1 Obtenció de les components principals

# Seleccionem les dades
X<-wine[,2:14]
# Centrem i escalem cada columna
X<-apply(X,2,function(x) (x-mean(x))/sd(x)) # X<-scale(X) equivalentment
# Construïm la matriu A
A<-t(X)%*%X

Diagonalitzem la matriu i obtenim valors i vectors propis. Els coeficients de las components principals són els vectors propis.

P<-eigen(A)$vectors
rownames(P)<-colnames(wine)[2:14]
colnames(P)<-paste('PC',1:13,sep='')
P

##                      PC1          PC2         PC3         PC4         PC5
## Alcohol     -0.144329395 -0.483651548 -0.20738262 -0.01785630  0.26566365
## Malic_ac     0.245187580 -0.224930935  0.08901289  0.53689028 -0.03521363
## Ash          0.002051061 -0.316068814  0.62622390 -0.21417556  0.14302547
## Alcal_ash    0.239320405  0.010590502  0.61208035  0.06085941 -0.06610294
## Magnesium   -0.141992042 -0.299634003  0.13075693 -0.35179658 -0.72704851
## phenols     -0.394660845 -0.065039512  0.14617896  0.19806835  0.14931841
## Flavanoids  -0.422934297  0.003359812  0.15068190  0.15229479  0.10902584
## Nonflav_phe  0.298533103 -0.028779488  0.17036816 -0.20330102  0.50070298
## Proanthocy  -0.313429488 -0.039301722  0.14945431  0.39905653 -0.13685982
## Color_int    0.088616705 -0.529995672 -0.13730621  0.06592568  0.07643678
## Hue         -0.296714564  0.279235148  0.08522192 -0.42777141  0.17361452
## OD280.OD315 -0.376167411  0.164496193  0.16600459  0.18412074  0.10116099
## Proline     -0.286752227 -0.364902832 -0.12674592 -0.23207086  0.15786880
##                     PC6         PC7         PC8         PC9        PC10
## Alcohol      0.21353865  0.05639636  0.39613926 -0.50861912  0.21160473
## Malic_ac     0.53681385 -0.42052391  0.06582674  0.07528304 -0.30907994
## Ash          0.15447466  0.14917061 -0.17026002  0.30769445 -0.02712539
## Alcal_ash   -0.10082451  0.28696914  0.42797018 -0.20044931  0.05279942
## Magnesium    0.03814394 -0.32288330 -0.15636143 -0.27140257  0.06787022
## phenols     -0.08412230  0.02792498 -0.40593409 -0.28603452 -0.32013135
## Flavanoids  -0.01892002  0.06068521 -0.18724536 -0.04957849 -0.16315051
## Nonflav_phe -0.25859401 -0.59544729 -0.23328465 -0.19550132  0.21553507
## Proanthocy  -0.53379539 -0.37213935  0.36822675  0.20914487  0.13418390
## Color_int   -0.41864414  0.22771214 -0.03379692 -0.05621752 -0.29077518
## Hue          0.10598274 -0.23207564  0.43662362 -0.08582839 -0.52239889
## OD280.OD315  0.26585107  0.04476370 -0.07810789 -0.13722690  0.52370587
## Proline      0.11972557 -0.07680450  0.12002267  0.57578611  0.16211600
##                    PC11        PC12        PC13
## Alcohol     -0.22591696 -0.26628645 -0.01496997
## Malic_ac     0.07648554  0.12169604 -0.02596375
## Ash         -0.49869142 -0.04962237  0.14121803
## Alcal_ash    0.47931378 -0.05574287 -0.09168285
## Magnesium    0.07128891  0.06222011 -0.05677422
## phenols      0.30434119 -0.30388245  0.46390791
## Flavanoids  -0.02569409 -0.04289883 -0.83225706
## Nonflav_phe  0.11689586  0.04235219 -0.11403985
## Proanthocy  -0.23736257 -0.09555303  0.11691707
## Color_int    0.03183880  0.60422163  0.01199280
## Hue         -0.04821201  0.25921400  0.08988884
## OD280.OD315  0.04642330  0.60095872  0.15671813
## Proline      0.53926983 -0.07940162 -0.01444734

La nova matriu de dades s’obté.

newdata<-X%*%P
head(newdata)

##            PC1        PC2        PC3        PC4        PC5        PC6
## [1,] -3.307421 -1.4394023 -0.1652728 -0.2150246 -0.6910933  0.2232504
## [2,] -2.203250  0.3324551 -2.0207571 -0.2905387  0.2569299  0.9245123
## [3,] -2.509661 -1.0282507  0.9800541  0.7228632  0.2503270 -0.5477310
## [4,] -3.746497 -2.7486184 -0.1756962  0.5663856  0.3109644 -0.1141091
## [5,] -1.006070 -0.8673840  2.0209873 -0.4086131 -0.2976180  0.4053761
## [6,] -3.041674 -2.1164309 -0.6276254 -0.5141870  0.6302409 -0.1230834
##              PC7         PC8         PC9        PC10       PC11          PC12
## [1,] -0.59474883 -0.06495586 -0.63963836  1.01808396 -0.4502932  0.5392891439
## [2,] -0.05362434 -1.02153432  0.30797798  0.15925214 -0.1422560  0.3871456499
## [3,] -0.42301218  0.34324787  1.17452129  0.11304198 -0.2858665  0.0005819316
## [4,]  0.38225899 -0.64178311 -0.05239662  0.23873915  0.7574476 -0.2413387757
## [5,] -0.44282531 -0.41552831 -0.32589984 -0.07814604 -0.5244656 -0.2160546934
## [6,] -0.40052393 -0.39378261  0.15171810 -0.10170891  0.4044444 -0.3783653606
##              PC13
## [1,]  0.066052305
## [2,] -0.003626273
## [3,] -0.021655423
## [4,]  0.368444194
## [5,]  0.079140320
## [6,] -0.144747017

Podem comprovar que les variables obtingudes estan incorrelacionades

round(cor(newdata),2)

##      PC1 PC2 PC3 PC4 PC5 PC6 PC7 PC8 PC9 PC10 PC11 PC12 PC13
## PC1    1   0   0   0   0   0   0   0   0    0    0    0    0
## PC2    0   1   0   0   0   0   0   0   0    0    0    0    0
## PC3    0   0   1   0   0   0   0   0   0    0    0    0    0
## PC4    0   0   0   1   0   0   0   0   0    0    0    0    0
## PC5    0   0   0   0   1   0   0   0   0    0    0    0    0
## PC6    0   0   0   0   0   1   0   0   0    0    0    0    0
## PC7    0   0   0   0   0   0   1   0   0    0    0    0    0
## PC8    0   0   0   0   0   0   0   1   0    0    0    0    0
## PC9    0   0   0   0   0   0   0   0   1    0    0    0    0
## PC10   0   0   0   0   0   0   0   0   0    1    0    0    0
## PC11   0   0   0   0   0   0   0   0   0    0    1    0    0
## PC12   0   0   0   0   0   0   0   0   0    0    0    1    0
## PC13   0   0   0   0   0   0   0   0   0    0    0    0    1

Les variàncies i desviacions típiques de cada component s’obtenen a partir dels valors propis

D<-eigen(A)$values
vars<-D/(nrow(X)-1)
round(vars,3)

##  [1] 4.706 2.497 1.446 0.919 0.853 0.642 0.551 0.348 0.289 0.251 0.226 0.169
## [13] 0.103

# La suma coincideix amb la suma de les variàncies de les dades originals
sum(apply(X,2,var))

## [1] 13

sum(vars)

## [1] 13

# El percentatge de variabilitat explicada per cada component s'obté
round(D/sum(D),3)

##  [1] 0.362 0.192 0.111 0.071 0.066 0.049 0.042 0.027 0.022 0.019 0.017 0.013
## [13] 0.008

Comprovem que s’obté el mateix a partir de les components principals

round(apply(newdata,2,var),3)

##   PC1   PC2   PC3   PC4   PC5   PC6   PC7   PC8   PC9  PC10  PC11  PC12  PC13 
## 4.706 2.497 1.446 0.919 0.853 0.642 0.551 0.348 0.289 0.251 0.226 0.169 0.103

1.15.2 Representació gràfica dels resultats

plot(D/sum(D),ylab='%',main='Percentatge de variabilitat aportat per cada
     component principal',cex.main=0.75)

# Representació dels individus segons les dues primeres components
plot(newdata[,1],newdata[,2],main='PCA wine data set',xlab='PC1',ylab='PC2',type='n')
text(newdata[,1],newdata[,2],label=wine[,1],cex=0.5,col=as.numeric(wine[,1]))
abline(h=0,col='red');abline(v=0,col='red')

Compareu els resultats amb el que obtindríem si fem les gràfiques de les observacions amb dues de les variables originals

plot(wine[,2],wine[,3],main='wine data set',xlab='Alcohol',ylab='Malic acid',type='n')
text(wine[,2],wine[,3],label=wine[,1],cex=0.5,,col=as.numeric(wine[,1]))

1.15.3 Biplot

# Representació de les variables
plot(P[,1],P[,2],main='PCA wine data set',xlab='PC1',ylab='PC2',type='n')
text(P[,1],P[,2],label=colnames(wine)[2:14],cex=0.75)
abline(h=0,col='blue');abline(v=0,col='blue')

1.15.4 Correlació de cada variable amb les diferents components principals

Amb les dades escalades la representació coincideix amb la del cercle de correlacions, cada variable representada pel coeficient de correlació amb les components principals dels eixos.

mcor<-cor(X,X%*%P)
round(mcor,2)

##               PC1   PC2   PC3   PC4   PC5   PC6   PC7   PC8   PC9  PC10  PC11
## Alcohol     -0.31 -0.76 -0.25 -0.02  0.25  0.17  0.04  0.23 -0.27  0.11 -0.11
## Malic_ac     0.53 -0.36  0.11  0.51 -0.03  0.43 -0.31  0.04  0.04 -0.15  0.04
## Ash          0.00 -0.50  0.75 -0.21  0.13  0.12  0.11 -0.10  0.17 -0.01 -0.24
## Alcal_ash    0.52  0.02  0.74  0.06 -0.06 -0.08  0.21  0.25 -0.11  0.03  0.23
## Magnesium   -0.31 -0.47  0.16 -0.34 -0.67  0.03 -0.24 -0.09 -0.15  0.03  0.03
## phenols     -0.86 -0.10  0.18  0.19  0.14 -0.07  0.02 -0.24 -0.15 -0.16  0.14
## Flavanoids  -0.92  0.01  0.18  0.15  0.10 -0.02  0.05 -0.11 -0.03 -0.08 -0.01
## Nonflav_phe  0.65 -0.05  0.20 -0.19  0.46 -0.21 -0.44 -0.14 -0.11  0.11  0.06
## Proanthocy  -0.68 -0.06  0.18  0.38 -0.13 -0.43 -0.28  0.22  0.11  0.07 -0.11
## Color_int    0.19 -0.84 -0.17  0.06  0.07 -0.34  0.17 -0.02 -0.03 -0.15  0.02
## Hue         -0.64  0.44  0.10 -0.41  0.16  0.08 -0.17  0.26 -0.05 -0.26 -0.02
## OD280.OD315 -0.82  0.26  0.20  0.18  0.09  0.21  0.03 -0.05 -0.07  0.26  0.02
## Proline     -0.62 -0.58 -0.15 -0.22  0.15  0.10 -0.06  0.07  0.31  0.08  0.26
##              PC12  PC13
## Alcohol     -0.11  0.00
## Malic_ac     0.05 -0.01
## Ash         -0.02  0.05
## Alcal_ash   -0.02 -0.03
## Magnesium    0.03 -0.02
## phenols     -0.12  0.15
## Flavanoids  -0.02 -0.27
## Nonflav_phe  0.02 -0.04
## Proanthocy  -0.04  0.04
## Color_int    0.25  0.00
## Hue          0.11  0.03
## OD280.OD315  0.25  0.05
## Proline     -0.03  0.00

plot(mcor[,1],mcor[,2],xlab='PC1',ylab='PC2',type='n',xlim=c(-1,1),ylim=c(-1,1))
text(mcor[,1],mcor[,2],labels=rownames(mcor),cex=0.75)
abline(h=0,col='blue');abline(v=0,col='blue')
nvar<-nrow(mcor)
arrows(rep(0,nvar),rep(0,nvar),mcor[,1],mcor[,2],col='red')